#include "analyzer_common.h"

void ilu0_ref(
    sparse_pointer_t* rowptr, 
    sparse_index_t* colidx, 
    sparse_value_t* lu, 
    sparse_index_t m, 
    sparse_pointer_t nnz, 
    sparse_index_t* diagptr
){
    for (sparse_index_t i = 0; i < m; ++i)
    {
        for (sparse_pointer_t j = rowptr[i]; j < diagptr[i]; ++j)
        {
        sparse_index_t c = colidx[j];
        sparse_value_t tmp = lu[j] /= lu[diagptr[c]];
        sparse_pointer_t k1 = j + 1, k2 = diagptr[c] + 1;
        while (k1 < rowptr[i + 1] && k2 < rowptr[c + 1])
        {
            if (colidx[k1] < colidx[k2])
            ++k1;
            else if (colidx[k1] > colidx[k2])
            ++k2;
            else
            {
            lu[k1] -= tmp * lu[k2];
            ++k1;
            ++k2;
            }
        }
        }
    } // for each row
}

// void bilu0_ref(
//     sparse_pointer_t* rowptr, 
//     sparse_index_t* colidx, 
//     sparse_value_t* lu, 
//     sparse_index_t m, 
//     sparse_pointer_t nnz, 
//     sparse_index_t* diagptr,
//     sparse_index_t block_order
// ){
//     for (sparse_index_t i = 0; i < m; ++i)
//     {
//         for (sparse_pointer_t j = rowptr[i]; j < diagptr[i]; ++j)
//         {
//             sparse_index_t c = colidx[j];
//             // double tmp = lu[j] /= lu[diagptr[c]];
//             dblas_3x3_gemm_a_eq_axb(
//                 &lu[j*block_order*block_order],
//                 &lu[diagptr[c]*block_order*block_order]
//             );
            
//             sparse_pointer_t k1 = j + 1, k2 = diagptr[c] + 1;
//             while (k1 < rowptr[i + 1] && k2 < rowptr[c + 1])
//             {
//                 if (colidx[k1] < colidx[k2])
//                 ++k1;
//                 else if (colidx[k1] > colidx[k2])
//                 ++k2;
//                 else
//                 {
//                     // lu[k1] -= tmp * lu[k2];
//                     dblas_3x3_gemm_c_minus_axb(
//                         &lu[j*block_order*block_order],
//                         &lu[k2*block_order*block_order],
//                         &lu[k1*block_order*block_order]
//                     );
//                     ++k1;
//                     ++k2;
//                 }
//             }
//         }
//         dblas_3x3_inv(&lu[diagptr[i] * block_order * block_order]);
//         for(sparse_pointer_t j = diagptr[i] + 1; j < rowptr[i+1]; j++){
//             dblas_3x3_gemm_b_eq_axb(
//                 &lu[diagptr[i] * block_order * block_order],
//                 &lu[j * block_order * block_order]
//             );
//         }
//     } // for each row
// }



void bilu0_ref(
    sparse_pointer_t* rowptr, 
    sparse_index_t* colidx, 
    sparse_value_t* lu, 
    sparse_index_t m, 
    sparse_pointer_t nnz, 
    sparse_index_t* diagptr,
    sparse_index_t block_order
){
    for (sparse_index_t i = 0; i < m; ++i)
    {
        for (sparse_pointer_t j = rowptr[i]; j < diagptr[i]; ++j)
        {
            sparse_index_t c = colidx[j];
            // double tmp = lu[j] /= lu[diagptr[c]];
            dblas_3x3_gemm_a_eq_axb(
                &lu[j*block_order*block_order],
                &lu[diagptr[c]*block_order*block_order]
            );

            double* a = &lu[j*block_order*block_order];
            register double a00 = a[0], a01 = a[1], a02 = a[2];
            register double a10 = a[3], a11 = a[4], a12 = a[5];
            register double a20 = a[6], a21 = a[7], a22 = a[8];
            
            sparse_pointer_t k1 = j + 1, k2 = diagptr[c] + 1;
            while (k1 < rowptr[i + 1] && k2 < rowptr[c + 1])
            {
                if (colidx[k1] < colidx[k2])
                ++k1;
                else if (colidx[k1] > colidx[k2])
                ++k2;
                else
                {

                    // register double* a = &lu[j*block_order*block_order];
                    register double* b = &lu[k2*block_order*block_order];
                    register double* c = &lu[k1*block_order*block_order];

                    register double b0, b1, b2;
                    register double c0, c1, c2;

                    b0 = b[0], b1 = b[3], b2 = b[6];
                    c0 = 0, c1 = 0, c2 = 0;
                    c0 += a00 * b0;
                    c1 += a10 * b0;
                    c2 += a20 * b0;
                    c0 += a01 * b1;
                    c1 += a11 * b1;
                    c2 += a21 * b1;
                    c0 += a02 * b2;
                    c1 += a12 * b2;
                    c2 += a22 * b2;
                    c[0] -= c0, c[3] -= c1, c[6] -= c2;

                    b0 = b[1], b1 = b[4], b2 = b[7];
                    c0 = 0, c1 = 0, c2 = 0;
                    c0 += a00 * b0;
                    c1 += a10 * b0;
                    c2 += a20 * b0;
                    c0 += a01 * b1;
                    c1 += a11 * b1;
                    c2 += a21 * b1;
                    c0 += a02 * b2;
                    c1 += a12 * b2;
                    c2 += a22 * b2;
                    c[1] -= c0, c[4] -= c1, c[7] -= c2;

                    b0 = b[2], b1 = b[5], b2 = b[8];
                    c0 = 0, c1 = 0, c2 = 0;
                    c0 += a00 * b0;
                    c1 += a10 * b0;
                    c2 += a20 * b0;
                    c0 += a01 * b1;
                    c1 += a11 * b1;
                    c2 += a21 * b1;
                    c0 += a02 * b2;
                    c1 += a12 * b2;
                    c2 += a22 * b2;
                    c[2] -= c0, c[5] -= c1, c[8] -= c2;
                    
                    ++k1;
                    ++k2;
                }
            }
        }
        dblas_3x3_inv(&lu[diagptr[i] * block_order * block_order]);
        for(sparse_pointer_t j = diagptr[i] + 1; j < rowptr[i+1]; j++){
            dblas_3x3_gemm_b_eq_axb(
                &lu[diagptr[i] * block_order * block_order],
                &lu[j * block_order * block_order]
            );
        }
    } // for each row
}